# high confidence clouds or high confidence cloud shadows or fill values
mask_medconf <- function(x){
# use base R function to convert numeric into bits
bs <- intToBits(x)
if ( ((bs[1]) | # cloud
# (bs[2]) | # low cloud confidence
(bs[3]) | # shadow
(bs[4]) | # snow
(bs[5]) | # water
(bs[6]) | # aerosol (only low quality)
(bs[8]) | # subzero
(bs[9]) | # saturation
(bs[10]) | # illumination
(bs[11]) | # saturation
(bs[13]) | # slope
(bs[14])) # vapor
== 1){
return("flag") } else {
return("valid")
}
}
# due to the nature of the intToBits function it is not possible to perform
# this step in a mutate pipe. Therefor a vector is created via lapply,
# joined to then joined to the fs_LND df and finally flaged pixels can be removed
mask_idex <- lapply(raw_fs$QAI, mask_medconf) %>% unlist()
fs_LND_filtered <- raw_fs %>%
mutate(QAI = mask_idex) %>%
filter(QAI == "valid")Quite important to note here is that while the absolute number of pixels additionally filtered with more strict conditions applied is quite small, increasing these conditions did yield effective improvements in removing additional outlyers in the data (particularly evident in direct pixel scatter plot comparisons between LND05/07/08). Therefor the deliberate inclusions/exclusion of certain filtering criterion could disproportionately sway the final results of sensor comparisons and feature spaces.
# Katja´s method to filter QAI values
raw_fs_KK <- raw_fs
raw_fs_KK$qa_cloud <- bitwAnd(bitwShiftR(raw_fs_KK$QAI, 1), 3)
raw_fs_KK$qa_shadow <- bitwAnd(bitwShiftR(raw_fs_KK$QAI, 3), 1)
raw_fs_KK$qa_snow <- bitwAnd(bitwShiftR(raw_fs_KK$QAI, 4), 1)
raw_fs_KK$qa_water <- bitwAnd(bitwShiftR(raw_fs_KK$QAI, 5), 1)
raw_fs_KK$qa_aerosol <- bitwAnd(bitwShiftR(raw_fs_KK$QAI, 6), 3)
raw_fs_KK$qa_subzero <- bitwAnd(bitwShiftR(raw_fs_KK$QAI, 8), 1)
raw_fs_KK$qa_saturation <- bitwAnd(bitwShiftR(raw_fs_KK$QAI, 9), 1)
raw_fs_KK$qa_zenith <- bitwAnd(bitwShiftR(raw_fs_KK$QAI, 10), 1)
raw_fs_KK$qa_illumination <- bitwAnd(bitwShiftR(raw_fs_KK$QAI, 11), 3)
raw_fs_KK$qa_slope <- bitwAnd(bitwShiftR(raw_fs_KK$QAI, 13), 1)
raw_fs_KK$qa_vapor <- bitwAnd(bitwShiftR(raw_fs_KK$QAI, 14), 1)
raw_fs_KK$quality <- 0
raw_fs_KK$quality[raw_fs_KK$qa_cloud != 0 |raw_fs_KK$qa_snow != 0 |raw_fs_KK$qa_shadow != 0] <- 1
raw_fs_KK <- raw_fs_KK[raw_fs_KK$quality == 0, ]print(paste("number of observations following SHS method for QAI filtering:", round(nrow(fs_LND_filtered)/6)) )## [1] "number of observations following SHS method for QAI filtering: 35109"
print(paste("number of observations following KK method QAI filtering:", round(nrow(raw_fs_KK)/6)) )## [1] "number of observations following KK method QAI filtering: 74693"
print(paste("The difference in included observation:", round(nrow(fs_LND_filtered)/nrow(raw_fs_KK), digits = 2), "% (n=", (nrow(fs_LND_filtered)-nrow(raw_fs_KK)), ")") )## [1] "The difference in included observation: 0.47 % (n= -237504 )"
As seen above the differences in the filtering methods and parameterization is very marginal, at ~1% difference in amount of included observations
##Additional data filtering and auxilirary varaible and indicity computation
fs_LND <- fs_LND_filtered %>%
# reduce file size for testing
#sample_n(10000) %>%
# filter out extreme values
filter(across(BLUE:SWIR2, ~ . > 0),
across(BLUE:SWIR2, ~ . < 10000)) %>%
mutate(date = as.Date(str_sub(scene, 1, 9),format = "%Y%m%d"),
month = lubridate::month(date),
year = lubridate::year(date),
month_year = paste0(month,"_",year),
doy = lubridate::yday(date),
sensor = str_sub(scene, -5, -1),
SWIR_ratio = SWIR2/SWIR1,
NDVI = ((NIR-RED)/(NIR+RED)),
NDTI = ((SWIR1-SWIR2)/(SWIR1+SWIR2))) # source: https://www.mdpi.com/2072-4292/10/10/1657/htm## Warning: Using `across()` in `filter()` is deprecated, use `if_any()` or
## `if_all()`.
## Warning: Using `across()` in `filter()` is deprecated, use `if_any()` or
## `if_all()`.
ggplot(fs_LND, aes(doy, NDVI, color=year, group=year)) +
geom_smooth() +
scale_colour_viridis_c(option = "D") +
theme_minimal()## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
Here it can be seen that over time, average NDVI maximum values tend to increase over time, particularly early in the year before the growing season. Over such a long time frame it is possible though that this chance could be attributed to more irrigation, climate change, or LULC rather been strictly being resultant of deviations in sensor spec?
ggplot(na.omit(fs_LND), aes(doy, NDTI, color=year, group=year)) +
geom_smooth() +
scale_colour_viridis_c(option = "D") +
theme_minimal()## `geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'
Similar story here, with NDTI decreasing over time. NDTI is also negatively associated with tillage, so the dip in the growing season is explained by more bare soil (tilled) locations being progressively covered by PV. Sudden decreases during the growing season indicate a tilling event, where as a slow decline at the end of the growing season are more likely to indicate PV senescence.
fs_LND_long <- fs_LND %>%
pivot_longer(cols=c("BLUE":"SWIR2"),names_to = "wavelength", values_to = "reflectance") %>%
as.data.frame() %>%
mutate(wavelength_num = case_when(wavelength == "BLUE" ~ 482,
wavelength == "GREEN" ~ 562,
wavelength == "RED" ~ 655,
wavelength == "NIR" ~ 865,
wavelength == "SWIR1" ~ 1610,
wavelength == "SWIR2" ~ 2200),
reflectance = reflectance/100) # scale to percentggplot(fs_LND_long, aes(wavelength_num, reflectance, color=year, group=year)) +
# geom_smooth(formula = y ~ s(x, bs = "cs", k=6)) +
stat_summary(fun=mean, geom="line", size = 1) + # draw a mean line in the data
scale_colour_viridis_c(option = "D") +
theme_minimal()When plotting a simple mean of spectra across years it can be seen that more recent years dont defaco have higher reflectance values across the entire electromagnetic spectrum
ggplot(fs_LND_long, aes(wavelength, reflectance, color=sensor)) +
# geom_jitter() +
geom_boxplot() +
scale_colour_viridis_d(option = "D") +
theme_minimal()Looks like lots of extreme values (even post QAI filter), but when looking at the next density plot section, you can that relatively speaking, these extreme outlying values are exceedinly rare and marginal.
# facet wrap by wavelength
ggplot(fs_LND_long, aes(reflectance, color=sensor, fill=sensor)) +
geom_density(alpha = 0.05) +
#geom_jitter() +
scale_colour_viridis_d(option = "D") +
scale_fill_viridis_d(option = "D") +
scale_x_continuous(expand = c(0, 0)) +
#scale_y_continuous(expand = c(0, 0), limits = c(0, 0.02)) +
theme_minimal() +
guides(col = guide_legend(nrow = 3))+
theme(legend.position = "bottom") +
facet_wrap(~wavelength)# facet wrap by sensor
ggplot(fs_LND_long, aes(reflectance, color=wavelength, color=wavelength)) +
geom_density(alpha = 0.05) +
#geom_jitter() +
scale_colour_viridis_d(option = "D") +
scale_fill_viridis_d(option = "D") +
scale_x_continuous(expand = c(0, 0)) +
#scale_y_continuous(expand = c(0, 0), limits = c(0, 0.02)) +
theme_minimal() +
guides(col = guide_legend(nrow = 3))+
theme(legend.position = "bottom") +
facet_wrap(~sensor)## Warning: Duplicated aesthetics after name standardisation: colour
ggplot(fs_LND_long, aes(x = reflectance, y = sensor, fill = stat(x))) +
geom_density_ridges_gradient(scale = 1.2, rel_min_height = 0.01) +
scale_fill_viridis_c(name = "reflectance", option = "D") +
scale_x_continuous(expand = c(0, 0), limits = c(0, 75)) +
theme_minimal() +
facet_wrap(~wavelength)## Picking joint bandwidth of 0.257
## Picking joint bandwidth of 0.276
## Picking joint bandwidth of 1.36
## Picking joint bandwidth of 0.416
## Picking joint bandwidth of 0.638
## Picking joint bandwidth of 0.51
ggplot(fs_LND_long, aes(x = reflectance, y = wavelength, fill = stat(x))) +
geom_density_ridges_gradient(scale = 3, rel_min_height = 0.01) +
scale_fill_viridis_c(name = "reflectance", option = "D") +
scale_x_continuous(expand = c(0, 0), limits = c(0, 75)) +
theme_minimal() +
facet_wrap(~sensor)## Picking joint bandwidth of 0.785
## Picking joint bandwidth of 0.397
## Picking joint bandwidth of 0.396
## Picking joint bandwidth of 0.456
## Picking joint bandwidth of 0.845
reference_spectra <- read.csv("data/feature_space/sli_gen_dark_soils_0p4.csv",
encoding = "UTF-8") %>%
mutate(SWIR_ratio = SWIR2/SWIR1,
NDVI = ((NIR-RED)/(NIR+RED)),
NDTI = ((SWIR1-SWIR2)/(SWIR1+SWIR2)))
ggplot(reference_spectra, aes(NDVI, SWIR_ratio, color = cover)) +
geom_point() +
scale_colour_viridis_d(option = "D") +
theme_minimal() # No dissagregation
ggplot(fs_LND, aes(NDVI, SWIR_ratio)) +
stat_bin_hex(bins = 300) +
geom_point(data=reference_spectra, aes(NDVI, SWIR_ratio), color ="red") +
scale_fill_continuous(type = "viridis") +
theme_minimal() +
scale_x_continuous(expand = c(0, 0), limits = c(-0.2, 1)) +
scale_y_continuous(expand = c(0, 0), limits = c(0, 1.2))## Warning: Removed 7 rows containing non-finite values (stat_binhex).
Here we can see that there are two dense zones in the spectral feature spaces. The smaller and fainter zone is predominately populated by landsat 7 pixels (clearly seen in the following plot)
ggplot(fs_LND, aes(NDVI, SWIR_ratio)) +
stat_bin_hex(bins = 200) +
geom_point(data=reference_spectra, aes(NDVI, SWIR_ratio), color ="red") +
scale_fill_continuous(type = "viridis") +
theme_minimal() +
scale_x_continuous(expand = c(0, 0), limits = c(-0.2, 1)) +
scale_y_continuous(expand = c(0, 0), limits = c(0, 1.2)) +
facet_wrap(~sensor)## Warning: Removed 7 rows containing non-finite values (stat_binhex).
Here we can see that the spectral library fits particularly well to Landsat 8 (probably do to the fact that is was created using lansat 8 pixels?).
ggplot(fs_LND, aes(NDVI, SWIR_ratio)) +
stat_bin_hex(bins = 100) +
geom_point(data=reference_spectra, aes(NDVI, SWIR_ratio), color ="red") +
scale_fill_continuous(type = "viridis") +
theme_minimal() +
scale_x_continuous(expand = c(0, 0), limits = c(-0.2, 1)) +
scale_y_continuous(expand = c(0, 0), limits = c(0, 1.2)) +
facet_wrap(~month)## Warning: Removed 7 rows containing non-finite values (stat_binhex).
Here we can also see that across all sensors, the spectral library seems to triangulate the data points better during growing season months.
Displayed with only point per reference class to improve interpretability
# FROM 1985-2000
fs_LND |> filter(year<2001) %>%
ggplot(., aes(NDVI, SWIR_ratio)) +
stat_bin_hex(bins = 100) +
geom_point(data=reference_spectra[c(1,3,25),], aes(NDVI, SWIR_ratio), color ="red") +
scale_fill_continuous(type = "viridis") +
theme_minimal() +
scale_x_continuous(expand = c(0, 0), limits = c(-0.2, 1)) +
scale_y_continuous(expand = c(0, 0), limits = c(0, 1.2)) +
facet_wrap(~year)# FROM 2001-2022
fs_LND |> filter(year>2000) %>%
ggplot(., aes(NDVI, SWIR_ratio)) +
stat_bin_hex(bins = 100) +
geom_point(data=reference_spectra[c(1,3,25),], aes(NDVI, SWIR_ratio), color ="red") +
scale_fill_continuous(type = "viridis") +
theme_minimal() +
scale_x_continuous(expand = c(0, 0), limits = c(-0.2, 1)) +
scale_y_continuous(expand = c(0, 0), limits = c(0, 1.2)) +
facet_wrap(~year)ggplot(fs_LND, aes(NDVI, NDTI)) +
geom_point() +
stat_bin_hex(bins = 250) +
geom_point(data=reference_spectra, aes(NDVI, NDTI), color ="red") +
scale_fill_continuous(type = "viridis") +
scale_x_continuous(expand = c(0, 0), limits = c(-0.1, 1)) +
scale_y_continuous(expand = c(0, 0), limits = c(-0.25, 0.8)) +
theme_minimal()# +## Warning: Removed 5 rows containing non-finite values (stat_binhex).
## Warning: Removed 5 rows containing missing values (geom_point).
# scale_x_continuous(expand = c(0, 0), limits = c(-0.2, 1)) +
# scale_y_continuous(expand = c(0, 0), limits = c(0, 1.2)) ggplot(fs_LND, aes(NDVI, NDTI)) +
stat_bin_hex(bins = 200) +
geom_point(data=reference_spectra, aes(NDVI, NDTI), color ="red") +
scale_fill_continuous(type = "viridis") +
theme_minimal() +
scale_x_continuous(expand = c(0, 0), limits = c(-0.2, 1)) +
scale_y_continuous(expand = c(0, 0), limits = c(-0.25, 1.2)) +
facet_wrap(~sensor)## Warning: Removed 3 rows containing non-finite values (stat_binhex).
ggplot(fs_LND, aes(NDVI, NDTI)) +
stat_bin_hex(bins = 100) +
geom_point(data=reference_spectra, aes(NDVI, NDTI), color ="red") +
scale_fill_continuous(type = "viridis") +
theme_minimal() +
scale_x_continuous(expand = c(0, 0), limits = c(-0.2, 1)) +
scale_y_continuous(expand = c(0, 0), limits = c(-0.25, 1.2)) +
facet_wrap(~month)## Warning: Removed 3 rows containing non-finite values (stat_binhex).
# old computationally exorbinant method of and plotting density plots
# # Get density of points in 2 dimensions.
# # @param x A numeric vector.
# # @param y A numeric vector.
# # @param n Create a square n by n grid to compute density.
# # @return The density within each square.
# get_density <- function(x, y, ...) {
# dens <- MASS::kde2d(x, y, ...)
# ix <- findInterval(x, dens$x)
# iy <- findInterval(y, dens$y)
# ii <- cbind(ix, iy)
# return(dens$z[ii])
# }
#
# fs_density <- as.data.frame(fs_LND) %>%
# filter(SWIR_ratio != "Inf") %>%
# mutate(density = get_density((.)$NDVI, (.)$SWIR_ratio, n = 500))
# ggplot(fs_density, aes(NDVI, SWIR_ratio, color = density)) +
# geom_point() +
# geom_point(data=reference_spectra, aes(NDVI, SWIR_ratio), color ="red") +
# scale_colour_viridis_c(option = "D") +
# theme_minimal() +
# scale_x_continuous(expand = c(0, 0), limits = c(-0.2, 1)) +
# scale_y_continuous(expand = c(0, 0), limits = c(0, 1.2)) +
# facet_wrap(~sensor)fs_LND |>
group_by( year, sensor) %>%
summarise(n_observation = n()) %>%
ggplot(., aes(year, n_observation, color=sensor)) +
geom_line() +
scale_fill_continuous(type = "viridis") +
scale_colour_viridis_d(option = "D") +
theme_minimal() ## `summarise()` has grouped output by 'year'. You can override using the
## `.groups` argument.
fs_LND_obs_overview <- fs_LND %>%
group_by(month, year, sensor) %>%
summarise(n_observation = n()) %>%
mutate(timestamp = zoo::as.yearmon(paste0(month,"_", year), format="%m_%Y"))## `summarise()` has grouped output by 'month', 'year'. You can override using the
## `.groups` argument.
ggplot(fs_LND_obs_overview, aes(timestamp, n_observation, color=sensor)) +
geom_line(size=1) +
#geom_smooth( ) +
scale_fill_continuous(type = "viridis") +
scale_colour_viridis_d(option = "D") +
scale_x_continuous(breaks=seq(1984,2022,2)) +
theme_minimal()The acute fluctuations in observation numbers seems to stem from the fact that during winter months there is a lower observation density (plotted below)
(same data with some different vis options)
ggplot(fs_LND_obs_overview, aes(x=month, y=n_observation, group=month,fill=sensor, color=sensor)) +
# geom_boxplot(alpha=0.3, notch = T) +
geom_violin(alpha=0.3) +
geom_jitter(alpha=0.3) +
#geom_density_ridges_gradient(scale = 1.2, rel_min_height = 0.01) +
scale_colour_viridis_d(option = "D") +
scale_fill_viridis_d(option = "D") +
theme_minimal() +
scale_x_continuous(breaks=seq(1,12,1)) +
scale_y_continuous(limits = c(0, 1000)) +
facet_wrap(~sensor)ggplot(fs_LND_obs_overview, aes(x=month, y=n_observation, group=month,fill=sensor, color=sensor)) +
geom_boxplot(alpha=0.3) +
geom_jitter(alpha=0.3) +
#geom_density_ridges_gradient(scale = 1.2, rel_min_height = 0.01) +
scale_colour_viridis_d(option = "D") +
scale_fill_viridis_d(option = "D") +
theme_minimal() +
scale_x_continuous(breaks=seq(1,12,1)) +
scale_y_continuous(limits = c(0, 1000)) +
facet_wrap(~sensor)ggplot(fs_LND_obs_overview, aes(x=month, y=n_observation, group=month,fill=sensor, color=n_observation)) +
geom_boxplot(alpha=0.5) +
# geom_violin(alpha=0.3) +
geom_jitter(alpha=0.9) +
#geom_density_ridges_gradient(scale = 1.2, rel_min_height = 0.01) +
scale_colour_viridis_c(option = "D") +
scale_fill_viridis_d(option = "D") +
theme_minimal() +
scale_x_continuous(breaks=seq(1,12,1)) +
scale_y_continuous(limits = c(0, 1000)) +
facet_wrap(~sensor)The density of observations in summer months are far greater (although more variable), While during winter months are seem to be constantly few observations. This trend is stable across sensors, so mostly likely can largely attributed to increased cloud/snow cover during winter months.
fs_LND_loose <- fs_LND_long %>%
# filter(sensor %in% c("LND07", "LND08")) %>%
mutate(obs_time_place = paste0(plyr::round_any(doy, 2, f = floor),"_", year,"_",POINT_ID)) %>%
group_by(obs_time_place) %>%
filter(n() > 6)Total number of observations that match the loose condition of being within ten days of each other: 2.8334^{4}
Days of year with available observations: 2, 3, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 30, 31, 36, 37, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221, 222, 223, 224, 225, 226, 227, 228, 229, 230, 231, 232, 233, 234, 235, 236, 237, 238, 239, 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255, 256, 257, 258, 259, 260, 261, 262, 263, 264, 265, 266, 267, 268, 269, 270, 271, 272, 273, 274, 275, 276, 277, 278, 279, 280, 281, 282, 283, 284, 285, 286, 287, 288, 289, 290, 291, 292, 293, 294, 295, 296, 297, 298, 299, 300, 301, 302, 303, 304, 305, 306, 307, 308, 309, 312, 313, 314, 315, 316, 317, 318, 319, 320, 321, 322, 323, 324, 325, 326, 327, 330, 331, 332, 333, 334, 335, 338, 339, 340, 341, 342, 343, 348, 349, 350, 351, 354, 355, 356, 357, 362, 363
for (i in 1:4) {
scens <- unique(fs_LND_loose$obs_time_place)[c((((i-1)*20)+1):((i)*20))]
fs_LND_loose_subset <- fs_LND_loose %>%
filter(obs_time_place %in% scens) |>
filter(sensor %in% c("LND07", "LND08"))
print(ggplot( fs_LND_loose_subset, aes(wavelength_num, reflectance, color=sensor, fill=sensor)) +
# geom_density(alpha = 0.05) +
geom_line() +
# scale_colour_viridis_d(option = "D") +
#scale_fill_viridis_d(option = "D") +
scale_x_continuous(expand = c(0, 0)) +
#scale_y_continuous(expand = c(0, 0), limits = c(0, 0.02)) +
theme_minimal() +
guides(col = guide_legend(nrow = 3))+
theme(legend.position = "bottom",
axis.text.y = element_text(angle = 45),
strip.text.y = element_text(size = 8, angle = 330)) +
facet_wrap( ~obs_time_place) )
}fs_LND_scatter <- fs_LND_loose |>
select(obs_time_place,sensor,wavelength,reflectance) |>
pivot_wider(names_from = sensor, values_from = c(reflectance)) |>
rename(reflectance_LND04 = LND04,
reflectance_LND05 = LND05,
reflectance_LND07 = LND07,
reflectance_LND08 = LND08,
reflectance_LND09= LND09)my.formula <- y ~ x
fs_LND_scatter_LND_7_8 <- fs_LND_scatter |> select(-reflectance_LND04,-reflectance_LND05, -reflectance_LND09) |>
na.omit()
ggplot(fs_LND_scatter_LND_7_8,aes(reflectance_LND07, reflectance_LND08)) +
stat_bin_hex(bins = 100) +
stat_smooth(method = "lm", formula = my.formula, color="red") +
stat_poly_eq(formula = my.formula,
aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
parse = TRUE,
label.x =1.2,
label.y =1.05) +
geom_abline(intercept=1,slope=1, linetype=2) +
scale_colour_viridis_d(option = "D") +
scale_fill_continuous(type = "viridis") +
theme_minimal() +
annotate("text", x=7, y=105,
label= paste0("R2 =",round((cor((fs_LND_scatter_LND_7_8$reflectance_LND07),
(fs_LND_scatter_LND_7_8$reflectance_LND08),
use="complete.obs")^2),2)), size=4.5, hjust=0) +
annotate("text", x= 7, y=100, hjust=0, size=4.5,
label= paste0("RMSE ==",round(sqrt(mean(((fs_LND_scatter_LND_7_8$reflectance_LND07)-
(fs_LND_scatter_LND_7_8$reflectance_LND08))^2,
na.rm=TRUE)),2)), parse=TRUE)+
annotate("text", x=7, y= 95, size=4.5, hjust=0,
label=paste0("Bias = ", round((mean((fs_LND_scatter_LND_7_8$reflectance_LND07), na.rm=TRUE) -
mean((fs_LND_scatter_LND_7_8$reflectance_LND08),na.rm=TRUE)),2)))+
annotate("text", x=7, y=90,size=4.5,hjust=0,
label=paste0("MAE = ", round(mean(abs((fs_LND_scatter_LND_7_8$reflectance_LND07)-
(fs_LND_scatter_LND_7_8$reflectance_LND08))),2)))ggplot(fs_LND_scatter, aes(reflectance_LND07, reflectance_LND08)) +
stat_bin_hex(bins = 100) +
stat_smooth(method = "lm", formula = my.formula, color="red") +
stat_poly_eq(formula = my.formula,
aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
parse = TRUE,
label.x =1.2,
label.y =1.05) +
geom_abline(intercept=1,slope=1, linetype=2) +
scale_colour_viridis_d(option = "D") +
scale_fill_continuous(type = "viridis") +
theme_minimal() +
facet_wrap(~wavelength,
scales="free"
) ## Warning: Removed 45534 rows containing non-finite values (stat_binhex).
## Warning: Removed 45534 rows containing non-finite values (stat_smooth).
## Warning: Removed 45534 rows containing non-finite values (stat_poly_eq).
kable(
fs_LND_scatter |>
select(-reflectance_LND04,-reflectance_LND05,-reflectance_LND09) |>
group_by(wavelength) |>
summarise(R2 = round((cor((reflectance_LND07),
(reflectance_LND08),
use="complete.obs")^2),2),
RMSE = round(sqrt(mean(((reflectance_LND07)-
(reflectance_LND08))^2, na.rm=TRUE)),2),
Bias = round((mean((reflectance_LND07), na.rm=TRUE) -
mean((reflectance_LND08), na.rm=TRUE)),2),
MAE = round(mean(abs((reflectance_LND07)-
(reflectance_LND08)), na.rm=TRUE),2))
) %>%
kable_styling(full_width = T)| wavelength | R2 | RMSE | Bias | MAE |
|---|---|---|---|---|
| BLUE | 0.04 | 4.35 | 0.01 | 1.97 |
| GREEN | 0.06 | 4.02 | 0.06 | 1.85 |
| NIR | 0.67 | 6.13 | -1.23 | 3.80 |
| RED | 0.20 | 3.97 | 0.23 | 1.87 |
| SWIR1 | 0.45 | 4.23 | 0.74 | 2.61 |
| SWIR2 | 0.51 | 3.32 | -1.21 | 1.94 |
fs_LND_scatter_LND_5_7 <- fs_LND_scatter |> select(-reflectance_LND04,-reflectance_LND08, -reflectance_LND09) |>
na.omit()
ggplot(fs_LND_scatter_LND_5_7, aes(reflectance_LND07, reflectance_LND05)) +
stat_bin_hex(bins = 100) +
geom_abline(intercept=1,slope=1, linetype=2) +
scale_colour_viridis_d(option = "D") +
scale_fill_continuous(type = "viridis") +
theme_minimal() +
stat_smooth(method = "lm", formula = my.formula, color="red") +
stat_poly_eq(formula = my.formula,
aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
parse = TRUE,
label.x =1.2,
label.y =1.05) +
annotate("text", x=7, y=105,
label= paste0("R2 =",round((cor((fs_LND_scatter_LND_5_7$reflectance_LND07),
(fs_LND_scatter_LND_5_7$reflectance_LND05),
use="complete.obs")^2),2)), size=4.5, hjust=0) +
annotate("text", x= 7, y=100, hjust=0, size=4.5,
label= paste0("RMSE ==",round(sqrt(mean(((fs_LND_scatter_LND_5_7$reflectance_LND07)-
(fs_LND_scatter_LND_5_7$reflectance_LND05))^2,
na.rm=TRUE)),2)), parse=TRUE)+
annotate("text", x=7, y= 95, size=4.5, hjust=0,
label=paste0("Bias = ", round((mean((fs_LND_scatter_LND_5_7$reflectance_LND07), na.rm=TRUE) -
mean((fs_LND_scatter_LND_5_7$reflectance_LND05),na.rm=TRUE)),2)))+
annotate("text", x=7, y=90,size=4.5,hjust=0,
label=paste0("MAE = ", round(mean(abs((fs_LND_scatter_LND_5_7$reflectance_LND07)-
(fs_LND_scatter_LND_5_7$reflectance_LND05))),2))) ggplot(fs_LND_scatter, aes(reflectance_LND07, reflectance_LND05)) +
stat_bin_hex(bins = 100) +
stat_smooth(method = "lm", formula = my.formula, color="red") +
stat_poly_eq(formula = my.formula,
aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
parse = TRUE,
label.x =1.2,
label.y =1.05) +
geom_abline(intercept=1,slope=1, linetype=2) +
scale_colour_viridis_d(option = "D") +
scale_fill_continuous(type = "viridis") +
theme_minimal() +
facet_wrap(~wavelength,
scales="free"
) ## Warning: Removed 49668 rows containing non-finite values (stat_binhex).
## Warning: Removed 49668 rows containing non-finite values (stat_smooth).
## Warning: Removed 49668 rows containing non-finite values (stat_poly_eq).
kable(
fs_LND_scatter |>
select(-reflectance_LND04,-reflectance_LND08,-reflectance_LND09) |>
group_by(wavelength) |>
summarise(R2 = round((cor((reflectance_LND07),
(reflectance_LND05),
use="complete.obs")^2),2),
RMSE = round(sqrt(mean(((reflectance_LND07)-
(reflectance_LND05))^2, na.rm=TRUE)),2),
Bias = round((mean((reflectance_LND07), na.rm=TRUE) -
mean((reflectance_LND05), na.rm=TRUE)),2),
MAE = round(mean(abs((reflectance_LND07)-
(reflectance_LND05)), na.rm=TRUE),2))
) %>%
kable_styling(full_width = T)| wavelength | R2 | RMSE | Bias | MAE |
|---|---|---|---|---|
| BLUE | 0.05 | 3.65 | -0.64 | 1.91 |
| GREEN | 0.10 | 3.45 | -0.82 | 1.93 |
| NIR | 0.70 | 5.33 | -1.93 | 3.21 |
| RED | 0.28 | 3.25 | -0.29 | 1.74 |
| SWIR1 | 0.50 | 4.20 | -1.57 | 2.71 |
| SWIR2 | 0.54 | 2.80 | -0.01 | 1.50 |
Overall the correspondence between both Landsat 5/7 & 7/8 seems to be pretty good? R^2 values are around 90% and the MAE around a tick above 2%, which doesn’t seem to drastic? In both comparison cases LND07 seemed to have (very) slightly lower reluctance values. Notable is the occasional occurrence of extreme outlyers across all the visable wavelengths in the spectra. Given that polluted pixels such as snow and clouds should have been removed through the QAI filters I am unsure what could be the source of these observed deviations (especially since they are not consistent across wavelengths).
Some float from the ideal 1-to-1 line can be likely explained by the two day buffer in doy matching between pixels. The relative infrequence of these extreme outlyers could be explained by a major LULC event such as mowing or harvesting occurring exactly in this potential two day difference in observation date? Yet this does not then adress why the large deviance in reflectance would only be observed for blue,green,red pixels and not NIR+ ones…
link <- "/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/X0060_Y0040/1984-2021_001-365_HL_TSA_LNDLG_NDV_TSI.tif"
tsi <- rast(link)
plot(tsi)
tsi_df <- as.data.frame(tsi)for (i in c(1:#length(names(tsi))
10)) {
# Step 1: Call the pdf command to start the plot
png(file = paste0("/data/Dagobah/greengrass/schnesha/03_RBF_interpolation_LND_1982-2021/ts_gif/plot_", i, ".png")) #, # The directory you want to save the file in
# width = 12, # The width of the plot in inches
#height = 12) # The height of the plot in inches
plot(tsi[[i]])
# Step 3: Run dev.off() to create the file
dev.off()
}